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Abstract 

Enhancing and controlling chaotic advection or chaotic mixing within liquid droplets 
is crucial for a variety of applications including digital microfluidic devices which 
use microscopic "discrete" fluid volumes (droplets) as microreactors. In this work, 
we consider the Stokes flow of a translating spherical liquid droplet which we per- 
turb by imposing a time-periodic rigid-body rotation. Using the tools of dynamical 
systems, we have shown in previous work that the rotation not only leads to one or 
more three-dimensional chaotic mixing regions, in which mixing occurs through the 
stretching and folding of material lines, but also offers the possibility of controlling 
both the size and the location of chaotic mixing within the drop. Such a control was 
achieved through appropriate tuning of the amplitude and frequency of the rotation 
in order to use resonances between the natural frequencies of the system and those 
of the external forcing. In this paper, we study the influence of the orientation of the 
rotation axis on the chaotic mixing zones as a third parameter, as well as propose 
an experimental set up to implement the techniques discussed. 

Key words: Chaotic advection. Chaotic mixing. Resonances, Control, 
Microfluidics, Droplet, Stokes flow. 
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1 Introduction 



The concept of chaotic advection, also referred to as Lagrangian chaos, was 
introduced some twenty- years ago in order to enhance mixing in laminar flows, 
the latter being mostly two-dimensional, time dependent incompressible flows 
(see [1] for a historical development). 

Several works have shown that the presence of chaotic advection in a flow 
drastically increases the transport of passive particles. This jump in transport 
properties has been quantifled by the diffusion coeflicient (see, e.g., the works 
of [2|[3] : see also |4]l5l[6]l7] for experimental studies and (Il7ll8j for comprehensive 
reviews). Chaotic advection is generally obtained by adding a degree of free- 
dom to an incompressible two-dimensional flow. This degree of freedom can 
take the form of either time dependence, [4l[5f6l[7j or a third spatial dimension 

PUD]. 

Due to the fact that microfluidic systems are characterized by low Reynolds 
numbers, typical flows in such devices are laminar and turbulence is inexistent. 
Thus, one has to turn to a different strategy to achieve mixing, and chaotic 
advection is often the most efficient way to accomplish this. 

Microfluidics can either use continuous streams as in the case of microchannels 
or individual droplets in the so-called "digital microfluidic devices". In the 
latter, mixing of multiple reagents takes place within "discrete" fluid volumes 
(droplets), thus offering the possibility of using a multitude of droplets with 
each droplet playing the role of a microreactor pT][T2] . 

Mixing in microfluidic devices using chaotic advection has recently attracted 
much attention. While there are many passive strategies based on altering the 
channel geometry for flows in microchannels, the use of active techniques based 
on forcing (see, e.g., p!3lll4lll5m6lll7j ) has also proved to be efficient, especially 
at low Reynolds numbers [18]. The combination of both passive and active 
methods has been explored as well p!8lll9ll20l[2T] . 

Mixing inside a drop subjected to a forcing (at low Reynolds number) has 
been studied extensively in the literature [9llQII22ll23ll24ll25ll26j . It has also been 
demonstrated experimentally using periodic forcing [271128] . 

This article builds upon the work of [29] which investigated the effect of peri- 
odic forcing on a translating droplet and its chaotic fluid flow regions. 
While previous works [9lllOI30ll31j have shown the presence of chaotic advection 
in three-dimensional bounded steady flows, we concentrate here on unsteady 
flows and use the added unsteadiness to manipulate the obtained chaotic be- 
havior through resonances [26ll32ll33 j . 

The physical system and the dynamical system which describes it are outlined 
in Sec. [2| and the numerical results are described in Sec. [3| We flrst recall that 
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the chaotic mixing zone can be monitored in both location and size, which we 
show quahtatively by displaying Liouvillian sections. We also quantify the size 
of the mixing zone and study its variation as the amplitude, frequency and ori- 
entation of the rotation vary. In the last section, we propose an experimental 
device capable of implementing the dynamics studied in this paper. 



2 Dynamical system 

2. 1 The dynamical system 

Consider a spherical Newtonian droplet which is itself immersed in an incom- 
pressible Newtonian fluid. The drop undergoes a translating motion as well 
as a rigid body rotation, as in the work of [10]. Furthermore, the droplet is 
assumed to be spherical and thus the interfacial tension large enough. We also 
suppose a very small Reynolds number, i.e. i?e <C 1. Given these assumptions, 
it follows that Stokes flow is a reasonable approximation, with the internal and 
external flows satisfying the boundary conditions at the surface of the droplet, 
namely the continuity of velocity and the tangential stress balance at the in- 
terface. In this paper, we further explore the addition of a degree of freedom 
to the above problem by making the amplitude of the rotation periodic in 
time as in [29j. We also assume that the (mean) amplitude, frequency and 
orientation of the rotation can be varied. 

The internal flow is a combination of a steady Hill vortex-type base flow and 
a perturbation which takes the form of an oscillating rigid body rotation. The 
location X of a passive tracer satisfles the dynamical system 

X = V (X, t) = Vo {X) + a^{t)u: x X, (1) 
where V denotes the velocity of the tracer. 

In Eq. ([T]), the base flow is denoted by Vq and the perturbation consists of a 
rotation having a time-periodic amplitude a^{t) and a oriented along the unit 
vector (b. We now select a moving Cartesian coordinate system translating 
with the center-of-mass velocity of the droplet. Let the unit vector point in 
the direction of the translation and the unit vector lie in the a> — e^— plane. 
The following non- autonomous dynamical system then follows: 

u = X = zx - a^{t)ujzy, (2) 

V = y = zy + a^{t) {u^x - cu^z) , (3) 

w = z = l-2x^ -2y^ - z^ + a^{t)ujxy, (4) 
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where all lengths and velocities have been made dimensionless by normalizing 
with respect to the droplet radius and the magnitude of the translational 
velocity. In this work, we assume that the amplitude of the rotation is purely 
sinusoidal about a mean value (it contains only one harmonic), i.e., 

ciuj{t) = - (1 + COS cut) with < 6: <C 1, (5) 



and that the orientation of the rotation is given by the unitary vector 

a> = {cjx^ 0, cJz) = (cos a, 0, sin a) . (6) 



Note that Eq. ([T]) is identical to the equations given in [lOj, except that the 
vorticity vector is no longer constant but now given by a^{t)Cj. This variation 
in time could result from the presence of unsteady vorticity in the external 
flow field or a time dependent body force. In practice, this could be realized, 
e.g., by creating a time dependent swirl motion in the external fiow or by 
applying an electric field capable of exerting a torque on the droplet (this 
could be achieved by using traveling wave dielectrophoresis which translates 
and rotates particles [34j or electrorotation which rotates particles [35]). Note 
that this fiow is the superposition of a Hill's vortex and an unsteady rigid 
body rotation, and that the surface of the droplet, = + + = 1, is 
invariant under the fiow given by Eqs. g, g and Q. 



2.2 Base flow 



In this section, we define the base fiow (^ = 0), which corresponds to the in- 
tegrable case in terms of dynamical systems theory. This is a two-dimensional 
axisymmetric fiow with two independent integrals of motion, e.g., the stream- 
function ifj and the azimuthal angle 0: 

= - [x^ + y^^ — r^^ , = arctan {y/x) ^ (7) 



where G [0, 1/8]. The streamlines, r^,(/,, are simply lines of constant and 0, 
satisfying the equation (1 — 2x^ — 2y^)^ + (2z(x^ + y^))^ = l — 8ip (see Fig. 2.2). 
The surface of the droplet coincides with heteroclinic orbits defined hjip = 
connecting two hyperbolic fixed points located at the poles of the spherical 
drop. As ip increases from ip = to the value ip = 1/8, the streamlines are 
closed curves converging toward a circle of degenerate elliptic fixed points 



(x^ + = 1/2, z = 0). The frequency of the dynamics on a streamline 
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Fig. 1. (a) Droplet in the perturbed case (with unsteady rigid body rotation); (b) 
Streamhnes within the droplet for the base (integrable) flow (without rigid body 
rotation). The motion on each streamline has a frequency Q (ip) given by Eq. Q. 

takes the expression 



27r 



7t/2 
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+ 7(^) sin/3 VI + 7 VVl + 7 



(8) 



where 7('0) = \/l — 8?/; and K refers to the complete elliptic function of the 
first kind. The frequency lies in between f](0) = and 1^(1/8) = \/2 (see 



Fig. 2.2). We now define a uniform phase x (27r) such that x = on 

the Bx — e^;— plane and X = ^ (ip). Every point in the interior of the droplet, 
except those lying on the z— axis, can be described by the values of (^0, (/), x). 
The unperturbed fiow can then be expressed in terms of the action-action- 
angle variables {ip^ 0, x) as follows: 

^ = 0, = 0, x = ^W- 



2. 3 Perturbed flow 



The perturbed fiow is defined by < ^ <C 1, with the integrals of motion ip 
and satisfying 

^ = -2a^{t)uJxip sin (pG{^p,x) , = a^{t)cj^ - a^{t)uJxCO^(pG {ip,x) , (9) 
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where x) = z/{x'^ + y'^) is a 27r periodic function of x, with zero average. 
The time evolution of x, on the other hand, satisfies the equation 



where i7 is a 27r periodic function of x- The perturbed system is characterized 
by two time scales, a fast one related to the evolution of x ^ slow one 
linked to both ip and cj). As we see below, chaotic advection can be obtained by 
exploring resonance phenomena between the frequency Q of the unperturbed 
flow (integral case) and the forcing frequency uu. 

From a dynamical systems viewpoint, the flow V{X^t) is the superimposi- 
tion of the integrable flow Vo{X) and a small time-dependent perturbation 
a^{t)u: X X, as given in Eq. ([T]). Since the unperturbed flow has two in- 
variants, the trajectories of this integrable dynamics are all periodic orbits. 
Most periodic orbits, however, are expected to break under the influence of 
a generic time-dependent perturbation with arbitrarily small amplitude, thus 
possibly leading to chaotic mixing properties. As well-known, the trajectories 
of an integrable system with only one invariant [36] are two-dimensional tori. 
Perturbing such an integrable system leads to poor mixing properties since 
two-dimensional tori (which are robust to small perturbations) act as barri- 
ers to chaotic diffusion. As in our previous work [29j, the goal of this paper 
is the generation of three-dimensional chaotic mixing regions of a given size 
and at speciflc locations. The strategy we adopt consists of bringing a family 
of unperturbed tori {F^^}^^^^ into resonance with the perturbation a^{t) by 
adjusting the frequency uj in order to satisfy the resonance condition: 

nVt{iljn) =cj. (10) 



The notation CMR^ is then used to denote the chaotic mixing region created 
around the torus F^^. 

We seek to control the mixing by varying the three parameters of the rotation, 
i.e., its amplitude 6:, its frequency u and its orientation a. Note that the effects 
of the amplitude and frequency have been studied for a flxed value of a [29j . 

The effect of a rotation is studied, with an amplitude such that 6: <C 1, with 
6: = at which chaotic mixing is inexistent and £ = e^nax at which mixing is 
maximum. In addition, we limit our study to the range < 6J < 72 which 
includes the frequencies ^ of all tori within the droplet. Notice that = on 
the boundary of the droplet and = \/2 as the torus approaches the elliptic 
fixed point. Furthermore, due to the symmetry of Eqs. ([2]), ([3]) and Q and 
without loss of generality, we restrict our study to0<a<7r/2. 



6 



3 Numerical results 



3. 1 Controlling the location of the chaotic mixing region 



In this section, we recall our previous findings on the effect of the amplitude 
and frequency of the rotation for the orientation a = 7r/4, as the infiuence of 
the orientation a on these results is studied below. 

Figures |3.1[ |3.1.1| and |3.1.1| display the Liouvillian sections of the perturbed 
flow, which consist of two-dimensional projections of time-periodic three-dimensional 
flows by a combination of a stroboscopic map and a plane section. Specifically, 
the Liouvillian sections considered here are the intersections of the trajectories 
with the plane ^ = at every period 2ti /uj. 

Figure 3.1 shows that the perturbation a^{t) creates two non-negligible three- 
dimensional chaotic mixing regions: one around the torus having the frequency 
uj denoted by CMRi and another one CMR^yi around the pole-to-pole axis 
and near the drop boundary. The latter region contains all tori having a fre- 
quency uj/n with n > 1. 

In Fig. 3J_, we clearly see that the location of CMRi varies with the value of 
uu according to Eq. (10). It is important to note that CMR^yi remains around 
the pole-to-pole axis and drop boundary due to the nearly vertical part of the 
curve close to ^0 = (see Fig. 2.2). 

For small values of cj, all resonances are located near the pole-to-pole hete- 
roclinic connections (at '0 = 0, near the z— axis and near the surface of the 
droplet, see Fig. 3.1). For larger uj values, CMRi separates from the chaotic 
region located close to the heteroclinic orbits and penetrates deeper into the 
droplet. In the interval < < \/2, CMRi is the largest chaotic region, 
followed by CMR^^i . As is increased further, CMRi moves toward the 
location of the central elliptic fixed point by following the location of the 
resonant torus with the frequency uj. 



3.1.1 Controlling the size of the chaotic mixing region 

In this section, we analyze the size of the two main chaotic mixing regions, 
i.e., CMRi and CMR^^i as the three parameters 6:, uj and a vary. In order to 
quantify the size, we use the fact that for a trajectory starting at ^0 = t/^q the 
adiabatic invariant varies between {ipo] a, uj) and {ipo] a, uj). It then 
follows that the width A?/; = ip^ {ipQ]e^a^uj) — {ip{)]e^a^uj) is large close 
to the resonance but decreases as one goes away from it. This quantity thus 
seems to be a good candidate to quantitatively estimate the size of CMRi 
and CMRnyi. 

As explained above, the location of CMRi is mostly determined by the rota- 
tion frequency uj (while CMR^^i is always located in the neighborhood of the 
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Fig. 2. Liouvillian sections for the rotation frequencies uj — 0.95, 1.1, 1.25, 1.40 (a-d), 
amplitude e — 0.05 and orientation a — 7r/4. The (red) hne inside the chaotic mixing 
region CMRi is the torus r^-i(^). 

heterochnic orbits). However, the size of CMRi and CMR^yi can be varied 
by adjusting the rotation amphtude £ and orientation a. This is illustrated 
in Fig. |3.1.1| which shows that the size of these chaotic mixing regions clearly 
increases with the amplitude of the perturbation. It is also interesting to note 
that around s = e^ax ^ 0.20, the two regions join and invade the entire drop. 
At that point, complete chaotic mixing is obtained. The size of CMRi as a 



function of the frequency uo is shown in Fig. 3.1.l| (middle panel). From this 
figure it is clear that for each value of e the size reaches a maximum for a 
certain value uj^ [e] of the forcing frequency. 

The effect of the orientation a on the size of CMRi is interesting. Indeed two 
distinguishable behaviors can be observed, one in which a hardly affects the 
size and another one in which a has a strong influence. These conclusions can 



be drawn from the lower panel of Fig. 3.1.1, where the size of the chaotic mix- 



ing region CMRi is displayed as a function of a. It is indeed clear that when 
a is not close to the limit values, i.e., and 7r/2, the size stays practically 
constant, but when a gets close to and 7r/2, it decreases significantly. Such 



observation is confirmed by looking at the Liouvillian sections in Fig. 3.1.1 
In the first two panels, we see a very small change in size while the size of the 
chaotic mixing regions starts decreasing in the third panel and shrinks drasti- 
cally in the fourth panel. Being aware of the dependence of the size of CMRi 



8 



1.0 
0.5 
0.0 
-0.5 

1.0 
0.5 
0.0 

Ay 

-0.5 

^ 

Fig. 3. (Color on line) Liouvillian sections for the rotation frequency u = 1.25, 
amplitudes e = 0.01,0.05,0.10,0.20 (a-d) and orientation a = 7r/4. The (red) line 
inside the chaotic mixing region CMRi is the torus r^-i(^). 

with respect to the orientation of the rotation a can be handy in practice in 
order to control the size of the mixing region. 

Indeed, in some applications where one is faced with the challenge of precisely 
controlling the size of the mixing zone despite uncontrollable, yet relatively 
small, fluctuations of a, setting a far away from the limit values and 7r/2 
and manipulating the chaotic zone size through the parameter £ should be 
desirable. In other, perhaps more gentle, applications where increasing the 
amplitude of the rotation may not be possible (e.g., in biomedical handling 
where biological particles need to be handled with care) , flxing s while tuning 
the size of the mixing by varying a around the limit values a = oi a = 7r/2 
could be the solution. One should notice, however, that at the critical orien- 
tation a = 7t/2 no chaotic mixing occurs and ip is conserved. 



4 Design of an experimental set up 

In this work, the flow within a drop was produced by the superimposition 
of an external steady translation and an unsteady rigid body rotation. The 
drop steady translating motion along the channel could simply be produced by 

9 




Fig. 4. Liouvillian sections for the rotation frequency = 1.25, amplitude e — 0.05 
and orientations a — 7r/64, 7r/8, 57r/16, Stt/S, Ttt/IG (a-f). The (red) hne inside the 
chaotic mixing region CMRi is the torus r^-i(^). 

means of a constant pressure gradient, exploring the buoyancy force in the case 
of a vertical column or using traveling wave dielectrophoresis in a microchan- 
nel ^J. The unsteady rigid rotation could be realized by applying an electric 
field that exerts a torque on the droplet as it is the case with traveling wave 
dielectrophoresis which generates both a force and a torque (although both are 
not independent) [34j. Another possibility is by using the so-called "electro- 
rotation" phenomenon which creates a torque [35j . Electrorotation stands for 
the spinning of an electrically polarized particle while the latter is subjected 
to a "rotating" electric field, that is an ac electric field generated by voltages 
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Fig. 5. Size of the chaotic mixing region. Upper panel: Normahzed size A'0 vs. 
orientation a for the amphtudes £ = 0.01,0.05,0.10,0.20 (a-d) and a frequency 
uj = 1.25; Middle panel: Normalized size Atp vs. frequency uj for the amplitudes 
e = 0.01, 0.05, 0.10, 0.20 (a-d) and an orientation o: = 7r/4 ; Lower panel: Normalized 
size A^/j vs. amplitude e for the frequencies u = 0.55,0.93,1.28,1.41 (a-d) and an 
orientation a = 7r/4. 



which are out of phase of one another. 

A possible design based on the latter phenomenon is proposed in Fig. [6) This 
device consists of a vertical square column/channel with electrodes embedded 
within its four walls creating a rotating electric field due to the phase difference 
between the voltages applied to adjacent electrodes in the same plane. It is 
known that such a four-pole electrode setting would generate a torque on the 
droplet trapped in the middle of the channel, with the torque strength being 
proportional to the square of the intensity of the electric field. The periodic- 
ity in the angular velocity is then obtained by imposing a phase diflFerence to 
the voltages applied between two four-pole electrode settings in two different 
X — ^—planes Pi and P^+i, so that a torque around the z— direction acts on 
the drop as it passes through Pi and a torque around the z— direction acts 
on the drop as it passes through P^+i. Although such a device would not pro- 
duce exactly the one harmonic angular velocity profile described in Eq. ([5]), 
our previous study [37] has shown that the strategy is robust with respect to 
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(a) (b) 




Fig. 6. Sketch of a possible design for the experimental apparatus showing (a) the 
square channel/column with a series of out-of-phase embedded electrodes, where 
uje and Ve stand for the frequency and amplitude of the voltage applied to the 
electrodes; (b) the four-pole electrode setting generating a torque in the z— direction, 
i.e., in the direction of translation, r > 0; (c) the four-pole electrode in a consecutive 
plane setting generating a torque in the z— direction, i.e. in the direction opposed 
to the direction of translation, r < 0. 

the specific time dependent function used for the rotation. Specifically, results 
were very similar when the sinusoidal rotation was replaced by a periodic tri- 
angular function. Finally, the control of the orientation of the rotation axis 
could be realized by a series of electrodes lying within inclined planes along 
the length of the channel (see Fig. [g]). 



5 Conclusions 



In this work, we have further studied the generation of chaotic mixing within a 
translating droplet by adding a perturbation in the form of an oscillatory rigid 
body rotation. As previously [29j, the frequency of the latter was selected in 
order to create resonances with the natural frequencies of the system, namely 
the frequencies of the various tori embedded within the drop. A particularly 
interesting feature of the perturbed system lies in the fact that both the size 
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and the location of the mixing region can be varied by adjusting the frequency 
and amphtude of the rigid body rotation. 

In this work, we have added a third parameter, namely the orientation of the 
rotation axis. It was found that the latter can influence the size of the chaotic 
mixing regions and that the size can signiflcantly decrease when the angle 
between the rotation axis and the direction of translation approaches either 
or 7r/2. Away from these two limit values, the size of the chaotic mixing region 
is maximal and the particular value of the angle has only a minor effect. 
A possible design for an experimental set up capable of guiding a drop in a 
controlled fashion by varying the amplitude, frequency and orientation of the 
rotation was also proposed. 
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